#ifndef __point_hh__
#define __point_hh__

#include <cassert>
#include <vector>

#include "femlib_exceptions.hh"

// Simpel class representing a n-D point
template <unsigned int DIM, class T=double>
class Point
{
public:
   T loc[(DIM == 0 ? 1 : DIM)];
   
   Point()
   {
      for (unsigned int i=0; i < DIM; ++i)
        loc[i] = 0;
   }
   
   Point(const T *ploc)
   {
      for (unsigned int i=0; i < DIM; ++i)
        loc[i] = ploc[i];
   }
   
   Point( const Point& pt )
   {
      for (unsigned int i=0; i < DIM; ++i)
        loc[i] = pt.loc[i];
   }

   template<typename U>
   Point( const Point<DIM,U>& pt )
   {
      for (unsigned int i=0; i < DIM; ++i)
        loc[i] = static_cast<T>(pt.loc[i]);
   }

   Point& operator=( const Point& pt )
   {
     if( this != &pt ) {
      for (unsigned int i=0; i < DIM; ++i)
        loc[i] = pt.loc[i];
      }
      return *this;
   }
   
   virtual ~Point()
   {
   }

   T& operator[]( unsigned int ind )
   {
     femlib_assert(ind < DIM, "Invalid index requested");
     return loc[ind];
   }

   T operator[]( unsigned int ind ) const
   {
     femlib_assert(ind < DIM, "Invalid index requested");
     return loc[ind];
   }

   void set(T *ploc)
   {
      for (unsigned int i=0; i < DIM; ++i)
        loc[i] = ploc[i];
   }
   
   void reset()
   {
      for (unsigned int i=0; i < DIM; ++i)
        loc[i] = 0.0;
   }
   
   T x() const
   {
     femlib_assert(DIM > 0, "DIM has to be greater than 0 to access 'x' value");
     return loc[0];
   }

   T y() const
   {
     femlib_assert(DIM > 1, "DIM has to be greater than 1 to access 'y' value");
     return loc[1];
   }

   T z() const
   {
     femlib_assert(DIM > 2, "DIM has to be greater than 2 to access 'z' value");
     return loc[2];
   }

};

// Simpel class representing a collection of n-D points
template <unsigned int DIM, class T=double>
class PointList
{
public:
   std::vector<T> loc;
   unsigned int n_points;
   
   PointList(unsigned int npts=DIM, T *ploc=NULL) : n_points(npts)
   {
     loc.resize(npts*DIM);
     if (ploc != NULL)
     {
      for (unsigned int j=0; j < DIM*n_points; ++j)
        loc[j] = ploc[j];
     }
   }
   
   ~PointList()
   {
     loc.clear();
   }

   PointList( const PointList& pt )
   {
      n_points = pt.n_points;
      loc.resize(n_points*DIM);
      for (unsigned int j=0; j < DIM*n_points; ++j)
        loc[j] = pt.loc[j];
   }

   PointList& operator=( const PointList& pt )
   {
     if( this != &pt ) {
      n_points = pt.n_points;
      loc.resize(n_points*DIM);
      for (unsigned int j=0; j < DIM*n_points; ++j)
        loc[j] = pt.loc[j];
      }
      return *this;
   }

   void set(T *ploc, unsigned int np)
   {
      femlib_assert (ploc != NULL, "The supplied pointer cannot be null") ;
      
      n_points = np ;
      loc.resize(n_points*DIM, 0);
      for (unsigned int i=0; i < DIM*n_points; ++i)
        loc[i] = ploc[i];
   }

   const T* get(unsigned int ipt) const
   {
     femlib_assert(ipt < n_points, "Invalid point index requested");
     return &loc[ipt*DIM] ;
   }
   
   T* get(unsigned int ipt)
   {
     femlib_assert(ipt < n_points, "Invalid point index requested");
     return &loc[ipt*DIM] ;
   }
 
   const Point<DIM,T> get_point(unsigned int ipt) const
   {
     femlib_assert(ipt < n_points, "Invalid point index requested");
     return Point<DIM,T>(&loc[ipt*DIM]) ;
   }

 
   Point<DIM,T> get_point(unsigned int ipt)
   {
     femlib_assert(ipt < n_points, "Invalid point index requested");
     return Point<DIM,T>(&loc[ipt*DIM]) ;
   }

   void reset()
   {
      for (unsigned int i=0; i < DIM*n_points; ++i)
        loc[i] = 0.0;
   }
   
   T x(unsigned int ipt) const
   {
     femlib_assert(DIM > 0, "DIM has to be greater than 0 to access 'x' value");
     return loc[ipt*DIM];
   }

   T& x(unsigned int ipt) 
   {
     femlib_assert(DIM > 0, "DIM has to be greater than 0 to access 'x' value");
     return loc[ipt*DIM];
   }


   T y(unsigned int ipt) const
   {
     femlib_assert(DIM > 1, "DIM has to be greater than 1 to access 'y' value");
     return loc[ipt*DIM+1];
   }
   
   T& y(unsigned int ipt) 
   {
     femlib_assert(DIM > 1, "DIM has to be greater than 1 to access 'y' value");
     return loc[ipt*DIM+1];
   }


   T z(unsigned int ipt) const
   {
     femlib_assert(DIM > 2, "DIM has to be greater than 2 to access 'z' value");
     return loc[ipt*DIM+2];
   }

   T& z(unsigned int ipt) 
   {
     femlib_assert(DIM > 2, "DIM has to be greater than 2 to access 'z' value");
     return loc[ipt*DIM+2];
   }


};

#endif
